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ABSTRACT 

We investigate the kinematic evolution of planetesimals in self-gravitating discs, com- 
bining Smoothed Particle Hydrodynamical (SPH) simulations of the disc gas with a 
gravitationally coupled population of test particle planetesimals. We find that at radii 
of 10s of au (which is where we expect planetesimals to be possibly formed in such 
discs) the planetesimals' eccentricities are rapidly pumped to values > 0.1 within the 
timescales for which the disc is in the self-gravitating regime. The high resulting ve- 
locity dispersion and the lack of planetesimal concentration in the spiral arms means 
that the collision timescale is very long and that the effect of those collisions that do 
occur is destructive rather than leading to further planetesimal growth. We also use 
the SPH simulations to calibrate Monte Carlo dynamical experiments: these can be 
used to evolve the system over long timescales and can be compared with analytical 
solutions of the diffusion equation in particle angular momentum space. We find that 
if planetesimals are only formed in a belt at large radius then there is significant scat- 
tering of objects to small radii; nevertheless the majority of planetesimals remain at 
large radii. If planetesimals indeed form at early evolutionary stages, when the disc 
is strongly self-gravitating, then the results of this study constrain their spatial and 
kinematic distribution at the end of the self-gravitating phase. 
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1 INTRODUCTION 

There are commonly discussed mechanisms by which gas 
giant planets can form out of the gas and dust of the pro- 
toplanetary disc. On the one hand, part of the disc may 
become sufficiently dense to become gravitationally unsta- 
ble and collapse (Boss 2000,). Alternatively, in what has be- 
come to be known as the 'core accretion' scenario, the dust 
in the disc coagulates to form ever larger objects, the plan- 
etesimals, eventually resulting in a solid core that is large 
enou gh to initiate r unawa y accretion of a gaseous envelope 
(e.g. iPollack et al.l l|l996h . iHubickvi et"aLri|2005l )'). In the 
latter scenario planetesimal coagulation must occur while 
the disc is still gas rich (in contrast to the case of the as- 
sembly of terrestrial planets, where the process may occur 
after dispersal of the disc gas). Unsurprisingly, therefore, 
there is a large literature devoted to p l anetes i mal evolution 
in the presence of gas (lAdachi et al.l (119761 ) IPoUack et al] 
l| 19961) ■ iTanaka fc WardI l|2004) '). bi most cases, the gas 
is treated as a laminar flow although several authors 
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(Lau ghlin eral] j2004) iNelson fc Papa loizou' ('2004l. lNelsoiil 
(2005f ). llda et al.[(|2008l ) j'Nelson fc GresseL (,2010 )) have con- 
sidered instead the scenario where planetesimals are sub- 
ject to stochastic torques arising from fluctuations in a disc 
subject to turbulence generated by the magneto-rotational 
instability (MRI). In the case of laminar discs, planetesi- 
mal eccentricity is (weakly) pum ped by mutual gravitational 
scattering (e.g. ISafronovl ([l96y)) and (weakly) damped by 
gas drag, so that the equilibrium distribution of eccentric- 
ities is peaked at low values (~ 0.01). This low value is 
important to continued planetesimal growth since it im- 
plies low relative velocities for planetesimal encounters: this 
not only f avours agglomerat i ve (a s opposed to dest r uctive ) 
outcomes (iBenz fc Asphaud (|l999l ). lLeinhardt et all (|2008l )) 
but also enhances the collision cross-section (oc e~^) in the 
gravitationally-focused regime. 

All of the above studies treat the evoluti on of plan- 
etesim als in non-self-gravitating discs. However, I Rice et aU 
(|2004 ) have argued that planetesimals may be formed very 
early (in the first few 10^ years) of a disc's life when it is 
still strongly self-gravitating; at this stage, spiral features in 
the disc provide pressure maxima in which dust is concen- 
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trated through the action of gas drag. The enhanced collision 
rates are favourable to grain growth and Rice et al. (2006) 
have argued that self-gravity in the dust phase may even 
promote the formation of km scale structures (i.e. planetesi- 
mals). Observational evidence for at least the early stages of 
grain growth during the self-gravitating phase is provided by 
the d etection of 10 cm radiation from HL Tau l|Greaves et al.l 
|2008| ). which implies that the growth of grains to at least cm 
scales (from the sub-micron scales typifying the interstellar 
medium) has already occurred in this young and massive 
disc system. 

If planetesimal formation indeed belongs to the earliest 
(self- gravitating) phase of disc evolution, then it is neces- 
sary also to trace the evolution of planetesimals during the 
self-gravitating phase. The dynamical evolution of planetes- 
imals in this environment has a number of implications for 
planet formation and for the coUisional production of dust 
in disc systems. For example, the relative importance of col- 
lisional growth of planetesimals in the self-gravitating and 
no n-self gravitati n g pha se of the disc can be crudely assessed 
fcf lBritsch erahl Hooi)) by comparing the product of the 
lifetime of each phase, the typical disc density and the (in- 
verse) square of the typical planetesimal velocity dispersion. 
The first two terms roughly cancel (i.e discs typically live an 
order of magnitude longer in the non-self gravitating phase 
but with disc masses an order of magnitude lower) so the 
relative importance of coUisional growth in the two regimes 
boils down to the relative values of the velocity dispersion. 

A pilot study of the dyna mical evolution of planetes- 
imals in self-gravitating discs l|Britsch et al.l [20081 ) demon- 
strated that high orbital eccentricities (e) are generated in 
such discs with particles undergoing stochastic changes in 
their orbital elements as a result of interaction with spiral 
features in the disc: high eccentricities imply a high veloc- 
ity dispersion (of order evk where vk is the local Keplerian 
velocity) if the particle trajectories are randomly phased. A 
lower velocity dispersion would however apply if the plan- 
etesimals instead demonstrated local v elocity coherence: the 
small number of particles modeled by iBritsch et al.l (|2008l ) 
did not permit exploration of this possibility however. For 
the same reason, it was not possible to discern the sign of 
any net migration of the planetesimal swarm nor whether 
the main statistical efi'ect involved changes in orbital energy 
or angular momentum. All these issues have implications for 
planetesimal growth during the self-gravitating phase and 
also for the reten tion of planetesima ls against loss through 
radial migration (Take uchi et al. * 2005) or coUisional grind- 
ing to small dust ( Wvatt et al.l 12007'). Moreover, evolution 
during the self-gravitating phase controls the spatial distri- 
bution and dynamical properties of any planetesimals that 
survive this phase. 

All of the above provides ample motivation for the 
present study in which we explore the response of a large 
ensemble of planetesimals to the gravitational torques im- 
posed by a self-gravitating gas disc. We restrict ourselves to 
the regime where planetesimals respond as test particles to 
the imposed potential fluctuations. This in practice restricts 
us to sizes > km s c ale (i n order to be able to neglect gas 
drag; iBritsch et all l|2008l ')) and < 1000 km scale (in order 
to be in the test particle regime where the efi'ect of gravita- 
tional perturba tions induced by t he particles in the disc gas 
can be ignored: [Bate et al.1 (120031 )'). 



In Setion 2 we present a numerical investigation (model- 
ing the disc with Smoothed Parti cle Hydrodynamics (S PH)) 
that extends the pilot study of IBritsch et al.l l|2008l ) to a 
large ensemble of planetesimals. Our results (Section 3) im- 
ply the disc pumps particle eccentricity with only mod- 
est changes in the semi-major axis and this motivates the 
analytical modeling and Monte Carlo simulations that we 
present in Section 4. This analytical/Monte Carlo approach 
is able to reproduce the results of the SPH simulation over 
the rather limited time frame that is possible in the lat- 
ter case and also permits integration over long timescales. 
In Section 5 we discuss the results of these calculations in 
relation to the questions posed above and summarise our 
conclusions. 



2 THE SIMULATION 

2.1 The physics of self-gravitating discs 

We model the gas disc as a self-gravitating disc which 
achieves a situation of self-regulation: i.e. it settles to a 
marginally-stable state where the Toomre Q parameter 
(Toomrc 1964) obeys 



ttGE 



(1) 



Here Cs is the local sound speed, E the gas surface den- 
sity and hi is the epicyclic frequency, which is ~ (the Kep- 
lerian frequency) for the low mass discs considered here. In 
this state the disc is in a state of thermal equilibrium be- 
tween the heating associated with the gravi tational instabil - 
ity and the imposed coolin g. Here we follow [Cammiel (|200ll ) 
and lLodato k, Ricd l|2005l ) in parameterising the cooling in 
terms of a cooling time that is a fixed multiple (/3) of the 
local dynamical time (Jl"^), i.e. we have fcooi = and 

Q-=-T^ (2) 

''cool 

where Q- and u are the cooling rate per unit mass and 
thermal energy per unit mass. 

The parameter /3 is a measure of overall cooling and the 
choice of j3 determines the magnitude of the spiral density 
waves required to maintain thermal equilibritum. By equat- 
ing this cooling law with the p redicted heating rate per unit 
mass Q+ from the instability, ICossins et al.l (|2009| ) demon- 
strated that the fractional density perturbation, a measure 
of relative spiral arm strength, is proportional to the 
numerical simulations in the same study confirmed this de- 
pendency and suggested that 
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(3) 



This means that in our constant /3 simulations we 
would expect the density perturbation to be approximately 
constant across the disc and to increase as /3 decreases. 
In practice, if /3 is less than about 4 in a simulation, 
the magnitude of t he perturbation becomes non-linear and 
the di sc fragments (lGanimiell200ll). (See e.g. iMeru fc Bate! 



( 201ll). iMerufcBatel (I2qi2l). iLodato fc Clarkd (l201ll) . 



Paardekooper et al. ( 201 j ) . jPaardekoopeJ { 20121 ) and 



Rice et alJ (|2012l )) for an ongomg debate as to whether frag- 
mentation may also eventually occur in well resolved simu- 
lations at significantly higher values of /3.) 
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2.2 Simulation parameters 

We model the system as a point mass, mass M, orbited 
by 250,000 SPH gas p articles and 50 , 000 te st particles; for 
details of the code see ICossins et all (|2009l ) and references 
therein. The test particles are assigned masses equal to 10~^ 
times the mass of a gas particle and are subject only to grav- 
itational forces. Gas particles are accreted onto the point 
mass if they enter within a si nk radius of 0.25 code units 
and satisfy certain conditions l|Bate et al.lll995l ): the point 
mass itself is free to move. By the end of the simulation, no 
more than 2 per cent of the gas particles are accreted. 

Artificial viscosity is included, according to the stan- 
dard SPH formalism, with qsph = 0.1 and /3sph = 0.2; 
for the parameters employed, the ratio of SPH smoothing 
length to disc vertical scale height is about two through- 
out. The gas is modelled as a perfect gas with 7 = 5/3 and 
is subject to both compressive heating and shock dissipa- 
tion. The cooling law (Equation [5)) means that heat lossq 
depends only on the dimensionless parameter j3, so cooling 
is scale free. For the main simulation we adopt P — 5 which 
corresponds to rather large amplitude spiral perturbations 
(see Equation [3]); we contrast these results with the lower 
amplitude case where /3 — 10. 

We employ scaled units so that G = M = 1 and ex- 
press time in units of the dynamical time tdyn = at 
R=l. This time t always refers to the time since the start of 
the simulation. We use these dimensionless results because 
they can be easily scaled to real situations. Our standard 
model consists of a disc with total mass equal to 0.1, dis- 
tributed with surface density oc R~^^^ over the radial range 
1 to 25 and the planetesimals are distributed according to 
the same radial distribution; we also investigate a similar 
set-up but with /3 = 10 and total mass of 0.5. We initialise 
the simulations with uniform temperature so that Q = 2 at 
all radii in the disc. At each radius we distribute the parti- 
cles in 2 in a Gaussian distribution representing approximate 
hydrostatic equilibrium at the initial temperature. True hy- 
drostatic equilibrium is reached within a few local dynamical 
times. 



3 SIMULATION RESULTS 

The disc is initialised with Q = 2 throughout and is thus ini- 
tially gravitationally stable. It then cools on the local cooling 
time until it attains Q 1, initiating the gravitational insta- 
bility and liberating heat. By time t — 4000 the disc has 
definitely settled into a quasis-steady state where spiral fea- 
tures form and dissolve on a roughly dynamical timescale 
and with Q 1 for 5 < R < 25. (At R < 5 the value of Q rises 
above unity since the disc is poorly resolved and heating 
by artifical viscosity acting on the Keplerian velocity field 
is sufficient to maintain the disc in a gravitationally stable 
state.) Since this is a purely numerical artefact, we restrict 
our attention to the regime R> 5. 

3.1 Escapees 

The first consideration is whether the test particles of the 
swarm gained enough energy from their interactions with 
the disc to escape the system entirely. We assess this by 



calculating the energies of the test particles (for this purpose 
we approximate the potential due to the disc by simply - for 
each particle - adding the enclosed disc mass to the effective 
mass of the central object and ignoring the contribution from 
external particles; although this is not strictly correct in a 
disc system, this does not significantly impair our ability to 
differentiate between bound and unbound particles). 

We find that (once the disc has settled into a steady 
gravoturbulent state) around 200 particles (out of a total 
of 10"") escape during the subsequent 4000 time units. If 
we (somewhat aribitrarily) equate code units for mass and 
radius with IMq and 1 an then this implies an e-folding 
timescale of 3 x 10^ years (i.e. of order the self-gravitating 
lifetime of the disc). This e-folding timescale is likely to be 
an under-estimate for two reasons: the low value of 13 em- 
ployed in the simulations (which implies much more vigorous 
spiral activity than would be expected given the re latively 
long c oolin g timescales in rea l istic p rotostellar discs; IClarkd 
(2009^) and lStamatellos et al.l l|2007h ) and the assumption of 
constant disc mass (whereas in reality the disc mass would 
decline over such timescales) . Q On this basis, we conclude 
that the majority of planetesimals would in reality remain 
bound. 



3.2 Ring spreading 

Once the system has settled into the equilibrium state (i.e. 
after t = 3 — 4000) we identify rings of particles over a 
restricted radial range and analyse their spatial evolution. 
Figure[T]illustrates the rather rapid spreading after five local 
orbital periods of a ring of initial width 1, initially located 
at i? = 10 at t = 4000. We find that the fractional change 
in the centroid of the distribution and the change in the 
width of the distribution (normalised to the initial radius) 
is independent of the radius of the ring selected as long as 
the comparisons are conducted after the same number of 
local dynamical timescales. This is as expected, given that 
(for the constant /3 that we impose) the fractional amplitude 
of the gravitational perturbations should be independent of 
radius. Note that although the mode of the particle radii 
has moved in slightly, the mean has moved out. Figure [2] 
illustrates the evolution of the fractional standard deviation 
of a ring at i? = 20 selected at t = 3500 and illustrates 
the continued spreading of the ring over many local orbital 
times. 

We find that the rings spread in response to stochas- 
tic interaction with the disc and also that there is a small 
outward shift in the particle centroid for each ring. We can 
understand this behaviour by considering the evolution of 
the energy distribution of particles in a selected ring (Fig- 
ure |3}: the finite width of the initial distribution just reflects 
the finite width of the ring. It is notable that there is lit- 
tle evolution in the distribution of particle energies over the 
time period of t = 4000 to t = 7500, whereas ring spreading 
is still significant over this period (Figure [2]). Although this 

^ We found an approximately 3 times higher escape rate in the 
case of our massive disc simulation (with disc mass five times 
higher than the standard case), which demonstrates that - as 
expected - the rate of test particle escape would drop as the disc 
mass dccHncs. 
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Radius 

Figure 1. Ring profile for i? = 10 after five local periods and 
with initial width of 1 at time t = 4000. 
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Figure 3. Evolution of the energy distribution of the R=10 ring, 
selected at t=4000. 
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Figure 2. Evolution of the ring standard deviation at i? = 20. 
The ring is selected at t = 3500, when the disc is in the quasi- 
steady state, and then followed for the rest of the simulation. 
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Figure 4. Eccentricity evolution of the standard simulation. 



spreading can be partially attributed to phase mixing (given 
the finite eccentricities of the particles selected), this is not 
the only process at work, as is illustrated by Figure [5] which 
shows the secular increase of the particle eccentricities. This 
suggests that as an ensemble, the particles are undergoing 
little evolution in energy but are diffusing in angular momen- 
tum. This implies a growth of particle eccentricity which we 
quantify below for the entire disc. (It should be stressed that 
the constancy of energy - and hence semi-major a xis - is a 
prope rty of the local ensemble; as was noted in .Britsch et al.l 
individual particles undergo stochastic interactions 
in which they undergo modest changes in energy). 

3.3 Eccentricity growth 

Figure [4] depicts the eccentricity distribution for the entire 
disc ensemble over the duration of the simulation and con- 
firms the evolution towards higher eccentricity, as would be 
expected if the spatial spreading largely reflects the growth 
of particle eccentricity at fixed semi-major axis. Note that 
in the period before the disc reaches the quasi-steady state 
(i.e. between t = and t = i — 4000) the eccentricity distri- 



bution is not being pumped by the recurrent spiral features. 
There does not appear to be a sharp transition between the 
two regimes though. 

In Figure [5] we compare the evolution of the mean ec- 
centricity between the standard simulations and the high /? 
case (in which /3 differs by a factor of 2). We see that the 
timescale to attain a given mean eccentricity is roughly dou- 
bled in the higher /3 case. This is consistent with a picture 
in which the particles are responding diffusively to gravi- 
tational interactions with the disc: the diffusion timescale 
scales with the inverse of the diffusion coefficient and thus 
with the inverse square of the amplitude of the perturba- 
tions. Hence (from Equation [3]) we expect the timescale for 
eccentricity growth to scale linearly with /?, as is consistent 
with Figure (5] We note that the mean eccentricity increases 
smoothly after t — 500 and that, as with the distribution 
as a whole, we cannot identify a sharp transition when the 
quasi-steady state is reached. This means that we can con- 
sider the distribution to be a reasonable proxy for how it 
would appear had the disc been in the quasi-steady state 
throughout, a fact we use in Section 4 to scale our analyti- 
cal model and Monte Carlo simulations. 
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Figure 5. Evolution of the eccentricity means for the two simu- 
lations. 



3.4 The velocity dispersion 

For a swarm of particles with randomly orientated orbits, 
the velocity dispersion is approximately equal to the prod- 
uct of the eccentri city and the local Keplerian velocity 
l|Britsch et alj [20081 ). The strong growth in particle eccen- 
tricity would thus raise the velocity dispersion of the plan- 
etesimals with consequences for particle growth that we dis- 
cuss below. However, this would be an over-estimate of the 
velocity dispersion if there was a degree of local velocity co- 
herence within the particle swarm. We investigate this in 
Figure |S] by comparing the evolution of the local velocity 
component (a) and its three orthogonal components (com- 
puted within a patch of disc containing ~ 100 particles) with 
the product of the instantaneous mean eccentricity and the 
local Keplerian velocity (cTp = e^k). The figure illustrates 
that eiik is indeed a good measure of the local velocity 
dispersion, as is expected in the case of randomly phased 
elliptical orbits. In addition, the ratio between the radial 
and azimuthal dispersions is maintained at about 3:2. One 
can show from the epicyclic approximation for a cylindrical 
potential that the exp ected ratio between the d ispersions 
should be cr^ = 0.5a| (jBinnev fc Tremainelliooil '). i.e. that 
= or 3(7^ ^ 2an. 

The growth of the velocity dispersion has profound 
consequences. First, a high velocity dispersion suppresses 
the role of gravitat ional focusing in facilitating collisions 
l|Britsch et al.|[2008l ). In the limit that the disc scale height 
exceeds the mean separation between planetesimals the 
timescale for physical collisions between planetesimals of 
mass M, radius R and density p is 



^"■"^ Epn(l + AGM/{Ra^)) ' ' ' 

where Sp is the planetesimal surface density and a is the 
planetesimal velocity dispersion. The a dependence means 
that collisions are highly disfavoured. At 30 au from a solar 
mass star and using e = 0.2, the velocity dispersion is around 
1 kms~^. Using the above formula for 100 km planetesi- 
mals, the collision timescale is of order 10^ years! Smaller 
planetesimals would give even longer timescales. Since plan- 
etesimal growth cannot proceed without physical collisions, 
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Figure 6. Evolution of the velocity dispersion (cr) and its com- 
ponents in the radial (crj;), azimuthal (cr^), and vertical {cr^) di- 
rections, compared with the prediction for randomly orientated 
orbits (cp). 

this rules out significant planetesimal growth during the self- 
gravitating phase. 

In addition - even when collisions do occur - the col- 
lisional velocity is correspondingly high. The velocity re- 
quired to shatter plane tesimals by collision is a few 100 ms~^ 
(|Leinhardt et al.ll2008l ). meaning that collisions would tend 
to result in disruption rather than in growth. We can then 
use the above collision probabilities to estimate a rough up- 
per limit on the amount of small dust (i.e. in the observable 
regime of < 1mm) that could be generated by disruptive 
collisions between planetesimals. On the optimistic assump- 
tion that every planetesimal-planetesimal collision results in 
all its mass being liberated as dust, the fraction of mass in 
planetesimals that can be liberated is simply the ratio of the 
self-gravitating lifetime to the collision time, i.e. ~ lO"**. We 
therefore conclude that the role of planetesimal collisions 
in re-supplying small dust is insignificant during the self- 
gravitating phase. Once the disc enters the non-self gravi- 
tating phase, and gas drag reduces the equilibrium eccentric- 
ity level to around 0.01 (Kokubo fc Ida 2000), the collision 
rate rises by around a factor of 100 and at this stage ei- 
ther coUisional growth and/or dust production may become 
important. 



3.5 Planetesimal concentration in the arms 

The above estimates for planetesimal collision frequencies 
neglect any possible concentration of particles in spiral fea- 
tures. For the values of /? employed in the simulations, the 
gas surface density varies by no more than a factor of two 
between the arms and the inter-arm regions (as expected, 
since we have modeled a regime where the disc is close to - 
but not at - the fragmentation boundary). If the planetes- 
imals simply followed the gas, such an enhancement would 
have only a minor effect on planetesimal collision probabil- 
ities, even in the unlikely event that planetesimals spent all 
their time in long-lived regions of density enhancement. In 
fact we find that - instead of being preferentially concen- 
trated in the arms - the planetesimal surface density varies 
by no more than ~ 20% around the orbit in the standard 
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simulation, i.e. with amplitude much less than that of the 
gas IT 



sPL 



iBritsch et al.l l|2008t ) suggested that there were possible 
hints that the concentration of particles in the arms was 
more efficient in the case of more massive discs. In their 
simulation with a disc mass of half that of the point mass 
they identified epochs at which some test particles would 
settle into orbits that co-rotated with the spiral modes and 
conserved a Jacobi constant. We have analysed our 'massive 
disc' simulation and however find that, as in the standard 
simulation, there is no discernible influence on the plan- 
etesimal distribution: there is no significant concentration 
of planetesimals in the arms. 

3.6 A truncated planetesimal distribution 

We finally consider the case that the planetesimal distribu- 
tion is truncated, i.e. that it does not initially extend within 
a certain inner radius, Rin, and study the extent to which 
planetesimals are scattered into regions within it^jn- This 
choice is motivated by the suggestion of IClarke Sz Lodatd 
l|2009l) that planetesimals are only likely to form in self- 
gravitating discs at large radius (be yond a few IPs of au) . 
We also note that the simulations of iGibbons et all (|2012l ) 
suggest that the mechanism is probably restricted to radii 
greater than 20 au. At such radii the cooling t ime i s rela- 
tively short (corresponding to /3 < 10; IClarkd l|200g| )) and 
the amplitude of spiral disturbances is large enough for rapid 
concentration of dust in spiral structures (i.e. on timescales 
shorter than the (roughly dynamical) lifetimes of spiral fea- 
tures). 

We do not need to conduct new simulations for this 
scenario but - since the planetesimals are non-interacting - 
simply tag planetesimals that are located at radius i? > 10 
at t = 5000 and follow the evolution of their density dis- 
tribution. From Figure [71 one can see that the distribution 
rapidly relaxes, with about 10 per cent of the particles mov- 
ing inward by a time of t=5500, although with no discernible 
further evolution over the time-frame of t = 5500 to 7500. 
The initial relaxation is consistent with the fact that the 
typical particle eccentricity is already ^ 0.1 at t — 5000 
and so just represents the fact that particles near the 'edge' 
visit the inner region even without further orbital evolution. 
We revisit the issue of planetesimal diffusion into an inner 
cleared region when we construct our analytic model in Sec- 
tion 4. 



4 A SIMPLE SCATTERING MODEL 

The simulation results encourage the exploration of a sce- 
nario in which the evolution of the planetesimal swarm is 
primarily driven by changes in angular momentum rather 
than energy. This prompted us to attempt to reproduce our 
results by modelling the evolution of the swarm as a process 



^ Note that the initial proc ess of concentratin g solid material into 
spiral features described bv lRice et al. - which is required 

for the formation of planetesimals in the self-gravitating phase - 
instead involves the effect of gas drag on small particles, and this, 
by definition, is ineffective in the test particle regime considered 
here. 
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Figure 7. Evolution of the planetesimals initially beyond R=10 
at time t=5000. 



of diffusion in angular momentum space. We start by ob- 
taining analytic expressions for the diffusion coefficients, set 
up the resulting diffusion equation and use Green's func- 
tions to solve for the evolution of the angular momentum 
distribution (and the associated eccentricity and radial dis- 
tributions) starting from an arbitrary radial profile and ini- 
tially circular orbits. The resulting solutions however involve 
the computation of large numbers of Legendre polynomials 
and are not of practical use until the system has evolved 
well away from the initial delta function distribution of ec- 
centricities. We nevertheless obtain solutions for the eccen- 
tricity distribution that are a reasonable representation of 
the SPH results at the end of the simulation. We can also 
compute the subsequent evolution of the system and derive 
expressions for the final (equilibrium) distributions over ra- 
dius and eccentricity (which correspond to a uniform distri- 
bution in angular momentum at every energy). In the fol- 
lowing sub-section we then verify these predictions using a 
simple Monte Carlo model; this approach has the additional 
advantage that it can be readily generalised to treat the case 
where the perturbation amplitude is a function of radius. 



4.1 Analytic solution 

We consider an ensemble of particles orbiting a mass M with 
fixed semi-major axis a. At some random time in each parti- 
cle's orbit the velocity vector is rotated by some angle A9 in 
the orbital plane, where is the angle between the velocity 
and radial vectors. This preserves the particle energy but 
perturbs the angular momentum. 

We start by writing the angular momentum per unit 
mass, L, for a particle at radius r as 



= 2GM sin^ e(r 



2a' 



(5) 



The equation is simplified by letting x = 1 — r /a, so 



that 

= GMa{l - x^) sm 
and 

di 



L2. 



(6) 



(7) 
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imax is the maximum local angular momentum and is 
a function of x (i.e. for particles at given r/a and energy 
- hence speed - this maximum occurs when the velocity is 
purely tangential). The angular momentum diffusivity, D, 
is given by AL'^/t where AL is the change in angular mo- 
mentum associated with deflection A6 and r is the time 
between such deflections: for now we assume both A9 and 
T to be constant around the orbit. We will also consider the 
case where is independent of a and where r scales with 
the orbital period (oc a^^^). This is motivated by the wish 
to compare with the SPH simulations which are conducted 
with constant cooling time to dynamical timescale ratio and 
in which, therefore, the fractional amplitude of gravitational 
disturbances is independent of radius (see equation (3)). 

We can relate AL to A9 via Equation [7| so that the 
value of D at given x is proportional to 4^, that is to say 



D oc {^f oc GMa{l 



(8) 



The expectation value of D over an entire orbit (at fixed 
L and a) is obtained by averaging this quantity over time, 
using the relationship between time interval and r for an 
elliptical orbit. Thus 



dr _ VGMa 
dt ^ r 



y/e^ - (1 - r/a)2 



(9) 



where e is the orbital eccentricity, which is related to L 
and a by 



= GMa{l ~ e^) 



(10) 



Equation [5] is derived from the orbital equation in the 
appendix. As before, the equation is simplified by subsitut- 
ing X for r, so that 



dt oc 



: da;. 



(11) 



(12) 



y/{^-x-) 
Thus we have that (for fixed a) 

where the odd terms integrate to zero over symmetric 
limits and may be discarded. The remaining even terms are 
standard integrals (requiring the substitution x — esinui) 
and yield 



{D) oc -(GMa 



(13) 



One can substitute this into Pick's Laws of diffusion 
to get the following ID diffusion equation for the angular 
momentum distribution h{L,t). J is the flux and A{a) is a 
constant that is determined by the strength of the pertur- 
bation and which in general depends on a (see below). In 
the case that we are currently considering (where the frac- 
tional amplitude of perturbations in the disc is independent 
of radius and where their characteristic timescale r scales 
with the orbital period i.e. oc a"^^^) then A is proportional 
to a~^^^; under these circumstances the timescale for angu- 
lar momentum diffusion scales with the local orbital period. 

Thus we have 



J : 



(14) 



dh 



dJ_ 



(D) = A{a){GMa- L'^). 



(15) 
(16) 



The form of (D) contrasts with other parametisations 
in the literature where it increases with L.|Ad ams fc BlochI 
(2009) consider, for example, the case where it is propor- 
tional to L. The form of Equation 1161 reflects the fact that 
where perturbations change only the direction of the veloc- 
ity, the associated angular momentum change goes to zero 
in the limit of tangential orbits; this form thus prevents dif- 
fusion of angular momentum beyond the physical limit set 
by a circular orbit. 

The diffusion equation is then 



dh 
'dt 



d_ 

dL 



A{a){GMa- 



dh 



for which we assume 
h{L,t) = M(L)T{t) to get 



(17) 



a separable solution: 



1 d 



{A{a){GMa~ 



dM, 



(18) 



1 dr 

T~dt A/ dL " ' dL 

This gives T{t) oc exp(—kt). If k is non-negative there 
will be no growing solutions. The equation for M{L) is 



dL 



((GMa - 



k 



Aia 



■M = 0. 



(19) 



Making the substitutions y — L/(\/ GMa) and 
k/A{a) = n(n + 1) converts Equation [19] into the Leg- 
endre equation: 



^((l-y')^) + n(n + l)M = 0. 



(20) 



The solutions must be regular at y = ±1, so we may 
use the Legendre polynomials. If this is combined with the 
solution for t we have 



h„(L, t) = c„exp(-n(n + l)A{a)t)P„{L/VGMa). 



(21) 



The general solution is a linear combination of the hn 
solutions, with the coefficients determined by the initial dis- 
tribution: 



h{L,t) = ^ Cnexf 



K;p(-n(n + l)A{a)t)PniL/VGMa). (22) 

n=0 

If the initial distribution consists of 
particles in prograde circular orbits, we have 
h{L,t = 0) = N5{L - VGMa). If we use the orthogo- 
nality property of the Legendre polynomials and the fact 
that -Pn(l) = 1 for all n, we have 



h{L,t) = 



N yr^2n + l 



VGMa 



exp{-n{n+l)At)P„{L/VGMa).{23) 



The eccentricity distribution is then obtained by using 
/(e)de = h{L)dL + h{~L)dL. This is because an eccentricity 
of e may be owing to a prograde or retrograde orbit with 
angular momenta of the same magnitude. The quantity dL 
is double- valued for the same reason. Thus from Eguation llOl 
we have: 



L = ±y/GMa{l-e^ 
and 



(24) 
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Figure 8. Evolution of the normalised eccentricity distribution. 



aL — — , de. 



(25) 



The Green's function solution, f{e,t,a) for initially cir- 
cular orbits at a particular semi-major axis a is therefore 



f{e,t,a) 



{2n+^)exp{-n{n+l)A{a)t)P„{^/^ - e^) 



To consider the behaviour of an entire disc of par- 
ticles we then calculate an appropriate superposition of 
Green's functions, bearing in mind (as argued above) that 
A = Aoa-"^''^. For comparison with the SPH simulations, 
we take a number distribution per unit a that is propor- 
tional to a~^/^ (since this corresponds to the initial surface 
density distribution scaling as R~^^^) and adopt inner and 
outer radii of ai = 1 and a2 — 25. The integrations are 
straightforward but laborious, making Mathematica the ob- 
vious choice. We plot F{e, t) (the distribution function for 
particle eccentricity) in Figure [8] and have normalised the 
time unit so as to match the peak in the distribution ob- 
tained from the SPH simulation at a time f=5000. At t=Q, 
F{e, t) is a delta function at e=0 but this and the subse- 
quent evolution before t=5000 required too many Legendre 
polynomials to be computationally practicable. We however 
note the similar form of the curve at t=5000 to that obtained 
in the SPH simulation at the same time (Figure 7). The fi- 
nal distribution occurs when t is large enough such that all 
the non-constant Legendre terms are effectively zero, leaving 
F{e,t) cx e/Vl — e^. This is an equilibrium state and repre- 
sents a uniform angular momentum distribution for particles 
with the same semi-major axis. 

The surface density distribution can be obtained 
through similar means. Equation [11] gives the probability 
of a particle with eccentricity e being at x. Multiplying this 
by /(e, t, a) and integrating with respect to e gives n{x, t, a), 
the fraction of particles at given semi-major axis in a given 
interval of x.The lower limit is the value of e for an orbit that 
achieves the required x value at either periapsis or apoapsis. 



n{x, t, a) oc 



f(e,t) 



Ve 



:de. 



(27) 



' 7nax(x , — x) 

This must be converted to a function of r and a and 
then normalised so that it represents the distribution re- 
sulting from a fixed number of particles. It can then be in- 
tegrated with the power law distribution for a to get the 
overall number density. 



N{r,t) (X / n{r,t,a)a da. 



S(r,t) = 



N{r,t) 
2Tvr 



(28) 



(29) 



Again, the integrations are complicated for very many 
Legendre polynomials. However the solution for the equilib- 
rium distribution is readily obtained. We note that (for a 
particle swarm of fixed a), the time spent in the interval r 
to r + dr is (^)^^dr, which (using Equation [9] and Equa- 
tion O can be written (for particles of angular momentum 
L) as 



dt oc 



rdr 



^rn — 

V -'^max 



L2 



(30) 



Integrating over all values of L from — Lmax to Lmax we 
obtain (in general): 



dt oc rdr 



h{L)dL 



(31) 



In equilibrium, h is constant which means that one may 
)te the integral in terms of the dimensionless variable L = 
L/Lmax SO that 



dt oc rdr 



dL 



Vl-L2 



(32) 



We thus find that in equilibrium the fraction of particles 
between r and r -I- dr scales as rdr. Now at fixed a, the 
minimum and maximum radii attained by the particle are 
and 2a so that the normalisation is such that 



T.^q(r,a) ^ < 2a) 



(33) 



We can now find the surface density density distribution 
by substituting n{r,t,a) = neq{r,a) into Equations 1281 and 
m\ Thus 



^, , A^(r,t) 1 
E(r,t) oc ^ ' ' oc - 
27rr r 



'^''^da(a > r/2). 



(34) 



This gives rise to three regimes. For < r < 2ai, all 
values of a contribute to the integral throughout this range 
and the surface density is constant. For 2ai < r < 2a2 then 
it is only a values greater than r/2 which contribute; conse- 
quently the surface density falls with radius as the annulus 
at r is populated by an ever decreasing range of a values. 
Finally, for r > 2a2, no particles can be scattered to this 
radius and the surface density is zero. 

Summarising we therefore have; 



E(r) = A(0 < r < 2ai), 

l-{ai/a2) 
E(r) = 0(r > 2r2). 



(35) 
(36) 
(37) 



We note that although we have derived the above for 
a specific distribution of particles in a, it is a general result 
that - in equilibrium - the surface density is uniform for radii 
less than twice the inner radius of the initial distribution. 
The reason for this is that every group of particles of given a 
is redistributed, in equilibrium, into a uniform density disc of 
particles extending from zero to 2a. At every radial position 
inward of 2ai , all bins of a provide uniform surface density 
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contributions to tiie total surface density. It is only at radii 
> 2ai that particles with a < r/2 cease to contribute and 
the surface density starts to decline. 

We emphasise that this entire derivation is based on the 
assumption that the disc fluctuations are scale-free, which 
implies that A9 is independent of a and independent of 
phase at given a. In reality, however, the amplitude of per- 
turbations in self-gravitating discs is an increasing function 
of radius. This changes the analysis in two ways: it changes 
the dependence of D on L (cf Equation I13p because the 
perturbations are of larger amplitude near apocentre and 
it also changes the scaling of D with a. As an example, if 
{d9)^ oc r^ then the diffusion coefficient becomes (by analogy 
with Equation I12p : 



(D) oc a" 



'' {GMajl ~ x^) - L^){1 - x) 



p+i 



■ dx. 



(38) 



In general, this integral is not analytic, except where 
p=-l: in this case the removal of the term (1 — x) makes 
no difference since the term oc x vanishes by symmetry for 
p = 0. Consequently the evolution at fixed a still reduces 
to the Legendre equation, although the computation for a 
range of a now needs to take account of the additional a 
dependence of D. We do not pursue this further here, since 
- in considering situations that are likely to occur in real 
discs - we want to be able to treat cases other than p — —1. 
This forms part of our motivation for adopting a Monte 
Carlo method in the following section. 



4.2 Monte Carlo simulations 

We have checked and extended the above analysis of the be- 
haviour of particle swarms that undergo diffusion in angu- 
lar momentum at fixed energy by conducting simple Monte 
Carlo simulations. Here we subject co-planar particles (ini- 
tially in circular orbits with the same surface density profile 
as employed in the SPH simulations) to random rotation of 
their velocity vectors through an angle in the range ±A6 at a 
random time each orbit. As expected, the evolution is quali- 
tatively independent of the value of A9 and with a timescale 
that scales as Ad~^. We have adjusted our parameters so as 
to match the peak of the eccentricity distribution in the SPH 
calculation at a time t — 5000 and verified that the evolu- 
tion of the eccentricity distribution (Figure [9]) is in good 
agreement with the output of our analysis of the diffusion 
equation (Figure [8)|. Secondly, the early evolution compares 
reasonably well with the f(e) curves obtained from the SPH 
simulation (Figure 3}. However, the equilibrium distribution 
for f(e) is only reached after several hundred thousand inner 
orbital periods. 

We also looked at the evolution of the surface density 
distribution of the planetesimals with initial radii greater 
that R=10 (Figure [To)) , so as to compare with that for the 
SPH simulation (Figure [7]). In the Monte Carlo simulation, 
the surface density profile evolves through a state similar to 
that seen in the SPH results and eventually attains a true 
equilibrium, where the surface density is approximately flat 
out to twice the inner radius (R=20)and declines to zero at 
twice the outer radius (R=50). This is the result expected 
from the earlier analysis. 
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Figure 9. Evolution of the un-normalised eccentricity distribu- 
tion from tlie Monte Carlo simulation. 
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Figure 10. Evolution of the surface density distribution, in ar- 
bitrary units, from the Monte Carlo simulation. 



4.3 Application to real discs 

We have hitherto adopted the assumption that the frac- 
tional amplitude of perturbations - and hence the value of 
A6 - is independent of radius. We have also normalised the 
evolution of the planetesimal swarm to that found in the 
SPH simulations which adopted rather short cooling times 
{P — 5 ~ 10) and correspondingly large amplitude of grav- 
itational instabi lity (with fractiona l amplitudes of several 
tens of percent jCossins et al.|[2010l )). This choice was mo- 
tivated by the need to compute the evolution of the plan- 
etesimal swarm in a reaonable time and also to ensure a 
sufficiently vigorous gravitational i nstability for it not to b e 
quenched by numerical viscosity ijLodato Sz Clark"3 l201lh . 
The insights provided by the simulations can now be used 
to rescale the problem to the parameter range expected in 
real self-gravitating discs. 

We therefore consider the case of an optically thick self- 
gravitating disc accreting at 3 x 1O~^M0 yr"'^ for which the 
steady st ate solutions (em ploying opacities due to dust and 
gas fromlBell fc LinI l| 1994 )1 are detailed in the Appendix of 
IClarkd l|2009l ). At radii > 35 au, the opacity is dominated 
by ice and we have: 
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tion from the Monte Carlo simulation with realistc cooling. 
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Figure 12. Evolution of the surface density distribution, in arbi- 
trary units, from the Monte Carlo simulation with realistc cooling. 



if ^ 



\ lOau 



-9/2 



(39) 



For 22au < R < 35au (where opacity is dominated by 
ice sublimation) and 



P = 270 



— y 

lOau J 



-9/14 



while for 4au < R< 22au 

-9/4 



/3 = 1000 



— 

lOau J 



(40) 



(41) 



Since the amplitude of pertubations scales as we 
also expect that A6 oc /3~^/^; we normalise the value of A9 
at /3 = 5 in order to match the rate of evolution of the 
SPH simulation. Within i? = 4 an we set AO = since it is 
arguable whether the disc is self-gravitating at this point; in 
any case, the long cooling times there equate to negligibly 
small values of A9. 

We start with a belt of planetesimals in circular or- 
bits which are distributed with a power law surface den- 
sity distribution with E oc since this is the steady 
state gas surface density profile for a self-gravitatin g disc 
with opacity dominated by ice grains (IClarkdliopgl '). The 
belt extends from 60 to 100 an, the choice of inner ra- 
dius being motivated by the fact that rather large ampli- 
tude density enhancements (rapid cooling) are required in 
order to for m planetesimals throug h dust concentration in 
spiral arms IClarke fc Lodatdl2009l i. A primary motivation 
for modeling a distribution with an initial hole is to discover 
whether - with a realistic prescription for the radial depen- 
dence of the perturbation amplitude - one expects planetes- 
imals to be scattered inwards to fill up the hole over the 
self-gravitating lifetime of the disc. 

We find (Figure I12|) that the evolution is qualitatively 
similar to that in Figure [TOl note that the edge is at 60 an in 
this simulation, as opposed to 10 code units in Figure [TOl In 
both cases, the time is normalised to the dynamical time at 
the inner edge (1 code unit and 1 au respectively); therefore 
in order to compare the two plots at a given number of or- 
bital times at the truncation point it is necessary to multiply 
the times in Figure [10] by 6^'^. It is then evident that the 



timescale for the infilling of the hole is rather similar in the 
two cases, i.e. the rate of infill is controlled by the ampli- 
tude of fiuctuations near the truncation radius. This can be 
understood inasmuch that particles that visit the inner disc 
on eccentric orbits spend most of their time near apocentre 
(i.e. close to the truncation radius). The fact that - in the 
realistic variable /3 case - the spiral structure is of very low 
amplitude in the inner disc therefore has little effect on the 
evolution of the planetesimal swarm. 



5 CONCLUSIONS 

The SPH simulations indicate that planetesimal eccentric- 
ity is efficiently amplified by interaction with spiral features 
in self-gravitating discs. In the case of discs where the cool- 
ing time is not much longer than the dynamical timescale 
(i.e in the outer disc, at radii > IDs of au, where we expect 
planetesimals to be formed in such discs), the amplitude of 
these features is sufficient for high eccentricities (> 0.1) to 
be driven on much less than the self-gravitating lifetime of 
the disc. We have found that there is no significant veloc- 
ity coherence within the particle swarm and thus that the 
local velocity dispersion is a significant fraction of the local 
orbital velocity; we also find that there is no tendency - in 
the test particle regime studied here - for the planetesimals 
to be concentrated within spiral arms. The lack of signifi- 
cant density enhancements and the large velocity dispersion 
(which weakens gravitational focusing) both contribute to 
very long collision times (> a Gyr); moreover, any collisions 
that did occur would be in the destructive regime. 

We have used the SPH simulations to calibrate Monte 
Carlo experiments in which the particle direction is ran- 
domly perturbed around its orbit and have compared these 
Monte Carlo experiments with an analytical description. 
The Monte Carlo experiments allow the long time integra- 
tion of the system and can also probe the high cooling time 
(weak spiral) regime that cannot be reliably simulated hy- 
drodynamically. We have treated the case of an initial plan- 
etesimal belt at large radius (60 au) where the radial varia- 
tion of the spiral amplitude is parametrised in terms of ex- 
pected variations in the local cooling physics in marginally 
unstable self-gravitating discs. We find that - notwithstand- 
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ing the fact that the spiral potential is very weak at small 
radius - there is significant scattering of planetesimals into 
the inner disc, with particles that are perturbed in the re- 
gion beyond 60 au visiting small radii on eccentric orbits. We 
nevertheless find that most of the particles initially beyond 
60 au are likely to be retained at large radius on timescales 
of ~ 10^ years. 

The picture that emerges from this study is that if 
planetesimals do form in the self-gravitating phase of disc 
evolution then they will be retained in the disc through- 
out this phase; their high eccentricities inhibit collisions 
and mean that they will neither grow nor suffer signifi- 
ca nt collisonal disru p tion over this period. If (as argued 
by IClarke fc Lodatd (|2009l )). such planetesimals are only 
formed in the outer disc then they will be largely retained 
at such radii, though with a significant minority that are 
perturbed to the inner regions of the disc. Such an endpoint 
would therefore represent the initial conditions for consid- 
ering the subsequent evolution of planetesimals during the 
non-self gravitating phase of disc evolution. 
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APPENDIX A: INTEGRATING AROUND AN 
ORBIT 

The expected value of D is obtained by averaging over one 
orbit. This requires obtaining dt as a function of r. We start 
with the orbital equation for r as a function of (j). 

a(l-e2) 



1 -I- e cos (f) 

dr dr d^ dr L 
dt ^ d^ dt ^ d0r2 

L is a constant of the motion for given values of a and 



(Al) 
(A2) 



L^ = GMa(l - e^ 



dr _ a(l - e^)esin0 v^GMa(l -e^) 
dt (1 -I- ecos</))2 r^ 

Substituting r back in simplifies the expression. 



(A3) 
(A4) 

(A5) 



dr _ VCMa sin 
dt ^ 1 - e2 a 

Re-arranging Equation I All and using that to eliminate 
(j) gives an expression in r and e only. 



dr _ %/GMa 
dt ^ r 



\/e2 - (1 - r/a)2 



(A6) 
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